molmass_acetate = 98; % Molar mass of potassium acetate
clc;
close all;

%% 37C acetate
tempmat = monod_mat_acetate_25C;
smat = size(tempmat);

ivec = [1:5]';

conc1 = tempmat(ivec,1)/100*1000/molmass_acetate*1000; % Convert to mM
grmax1 = nanmean(tempmat(ivec,2:4),2);
grmax_err1 = nanstd(tempmat(ivec,2:4),0,2);

figure;
errorbar(conc1,grmax1,grmax_err1, 'bo', 'linewidth', 1); 
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'b', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);

%% 30C acetate
tempmat = monod_mat_acetate_30C;
smat = size(tempmat);

ivec = [1:4]';

conc1 = tempmat(ivec,1)/100*1000/molmass_acetate*1000; % Convert to mM
grmax1 = nanmean(tempmat(ivec,2:4),2);
grmax_err1 = nanstd(tempmat(ivec,2:4),0,2);

hold on;
errorbar(conc1,grmax1,grmax_err1, 'ko', 'linewidth', 1); 
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'k', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);

%% 37C acetate
tempmat = monod_mat_acetate_37C;
smat = size(tempmat);

ivec = [1:4]';

conc1 = tempmat(ivec,1)/100*1000/molmass_acetate*1000; % Convert to mM
grmax1 = nanmean(tempmat(ivec,2:4),2);
grmax_err1 = nanstd(tempmat(ivec,2:4),0,2);

hold on;
errorbar(conc1,grmax1,grmax_err1, 'ro', 'linewidth', 1); 
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'r', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);